## reproduce the Nielsen analysis

## Required libraries: stm, gdata, igraph


## Set directories
## The directory with the original data
datadir <- "~/FinalPublicReplication/"
## A directory to output figures
outdir <- "~/FinalPublicReplication/figures/"
## The directory with the STM objects
indir <- "~/FinalPublicReplication/"



####################################################
## Run the model and label the topics
####################################################


## load the stm library
library(stm)

## set the working directory to where the STM output is saved
setwd(indir)
## load the Lucene object feeding into the STM model
load("CombinedLuceneREPLICATION.RData")
## run the model
#mod.out.15<-selectModel(data,vocab,K=15, prevalence=~combinedjihadi,data=meta, seed=02138, runs=20)
#save(mod.out.15,data,vocab,meta,file="RichCheck102314-SM.RData")
## load the STM object
load("RichCheck102314-SM.RData")

## This is the model we chose
mod.out <- mod.out.15$runout[[3]]
## save the topic number of the excommunication topic
eTopicNum <- 11
## save the topic number of the fighting topic
fTopicNum <- 1

## get the result of labelTopics
topicslists <- labelTopics(mod.out,n=10)

## My process of labeling topics was long.  I'm leaving the whole process here
## so that it's clear how I inferred the topic labels we use.

## First, I looked at the arabic results of labelTopics, focusing on the frequency and FREX and came up with some 
## quick labels that I thought might be right.
labelTopics(mod.out,n=10)

## Topics for mod.out.15$runout[[3]]
topicnames <- c("Fighting","Society","International","The Prophet","Prayer",
                "Ramadan","Family","Money, Pilgrimage, Marriage","Society???","Hadith",
                "Excommunication","Salafism","Shariah?","Cosmology","Hadith Narration")

## Then I started translating the Frequency and FREX terms carefully, both for use in the figure later and to 
## help label the topics.  There were some terms that I didn't know or wasn't sure about.  I marked these with
## question marks.

## Translation of frequency terms
englishF <- c("Muslim, Jihad, Islam, fight, Jihadi fighters, pathway, almighty, that, fighting, world", ## topic 1
              "person, life, soul/self, knowledge/science, society, work, PICTURE?, ?, thought, theory", ## topic 2
              "Arab, Jews, country, Islam, A.D., year, West, Muslim, religion, world", ## topic 3
              "said, prayers (be upon him), peace (be upon him), almighty, messenger, glory, prophet, that, almighty, people", ## topic 4
              "prayer, pray, son, prophet, sheikh, mosque, fatwas, group, theses", ## topic 5
              "day, fasting, Ashura, Ramadan, sheikh, group, fatwas, Uthaymeen, letters, theses", ## topic 6
              "woman, O???, man, girl, one, says, men, people, brother, day", ## topic 7
              "tithing, money, pilgrimage, permitted, religion, marriage, belief??, divorce, selling, forbidden", ## topic 8
              "Islam, land, mankind, people, religion, life, last/other??, al-???, ignorance, that", ## topic 9
              "Saying, hadith, said, prayers (be upon him), peace (be upon him), Muslim, legaly, not, thing/command?, prophet", ## 10
              "Apostasy, said, almighty, polytheism, Islam, Apostate, saying, people, faith, sheikh", ## 11
              "Sunna, sheikh, son, people, book, knowledge, Salafi, Muhammad, innovation, mercy", ## 12
              "Islam, wisdom, right, people, thing, legaly, Shari\'a, religion, descend, Muslim", ## 13
              "knowledge, characteristic, saying, meaning, ???, ???, people, said, faith, sn???", ## 14
              "Said, son, son, father, hadith, narrated, Abd, (may God be) pleased (with him), father, Omar") ## 15

## Translation of FREX terms
englishFREX <- c("jihad, fighting, jihadist fighters, pulpit, approves of us, ????, to fight, vicinity, martyr, factions", ## 1
                 "imagine, morals, develop, society, product, predestination?, environment?, traditions, activity, capitalism", ## 2
                 "capitol, Asia, Iran, South, Washington, A.D., Russia, Turkey, Mexico, Gulf", ## 3
                 "almighty, almighty, glory, bless you, witch?, punishment?, opinions?, sins?, fire, ???", ## 4
                 "prostration, prostrated, Abd al-Aziz (Bin Baz), supplicant, (Bin) Baz, prayer space, omission, prostration, mosque, Abdullah", ## 5
                 "wash, one who fasts, fasting, fasting, to break fast, Ramadan, travel, dirty, Ramadan retreat, ritual cleansing by dust", ## 6
                 "veil, youth, tamim??, Azzam?, tanks?, ??, ??, ??, sister, ??", ## 7
                 "tithing, divorce, banks, divorce, card, banks, to perform pilgrimage, poor, sell, ???", ## 8
                 "Europe, civilization, European, mankind, church, goods, generations, their lives, huge, legend", ## 9
                 "to forbid, analogy, permission, general, evidence, forbid, text, absolutely, voting, concensus", ## 10
                 "excommunicate, apostate, apostasy, ???, idolatry, excommunication, idols, to make permissible, ???, apostasy", ## 11
                 "heterodoxy, innovator, wool??, salafi, neighbors??, distinguish, ??, to go straight??, innovation, multiply", ## 12
                 "Shari\'a, to legislate, to send down, to judge, judgment, justice, parliament, court, democracy, those", ## 13
                 "creatures, characteristic??, names, throne, Tahawi??, ??, mean, ??, Ala\'?, explainer", ## 14
                 "narrate, Daoud, chain of narration, narrate, ???, Tabari??, our saying, Tarmidhi??, Father, Hurayra??") ## 15

## Then, to figure out the terms I wasn't sure about, I looked at the words in context.
## I made a search function that gets the top docs for each topic and searches them for an arabic term 
## (I used my own transliteration tools to make this easier).  This process relies on the raw
## texts which we are not releasing, so I cut out the code and just left the resulting translation.

## corrected frequency translations
englishF <- c("Muslim, Jihad, Islam, fight, Jihadi fighters, pathway, almighty, that, fighting, world", ## topic 1
              "person, life, soul/self, knowledge/science, society, work, image, material/physical, thought, theory", ## topic 2
              "Arab, Jews, country, Islam, A.D., year, West, Muslim, religion, world", ## topic 3
              "said, prayers (be upon him), peace (be upon him), almighty, messenger, glory, prophet, that, almighty, people", ## topic 4
              "prayer, pray, son, prophet, sheikh, mosque, fatwas, group, theses", ## topic 5
              "day, fasting, Ashura, Ramadan, sheikh, group, fatwas, Uthaymeen, letters, theses", ## topic 6
              "woman, O, man, girl, one, says, men, people, brother, day", ## topic 7
              "tithing, money, pilgrimage, permitted, religion, marriage, believe/ratify, divorce, selling, forbidden", ## topic 8
              "Islam, land, mankind, people, religion, life, other, God, ignorance, that", ## topic 9
              "Saying, hadith, said, prayers (be upon him), peace (be upon him), Muslim, legally, not, thing/command, prophet", ## 10
              "Apostasy, said, almighty, polytheism, Islam, Apostate, saying, people, faith, sheikh", ## 11
              "Sunna, sheikh, son, people, book, knowledge, Salafi, Muhammad, innovation, mercy", ## 12
              "Islam, wisdom, right, people, thing, legally, Shari\'a, religion, descend, Muslim", ## 13
              "knowledge, qualities, saying, meaning, Quran, to be, people, said, faith, Sunna", ## 14
              "Said, son, son, father, hadith, narrated, Abd, (may God be) pleased (with him), father, Omar") ## 15

## corrected FREX translation
englishFREX <- c("jihad, fighting, jihadist fighters, pulpit, approves of us, annotated, to fight, vicinity, martyr, factions", ## 1
                 "imagine, morals, develop, society, product,necessarily, environment, traditions, activity, capitalism", ## 2
                 "capitol, Asia, Iran, South, Washington, A.D., Russia, Turkey, Mexico, Gulf", ## 3
                 "almighty, almighty, glory, bless you, magic, punishment, hypocrisy, sins, fire, entreat", ## 4
                 "prostration, prostrated, Abd al-Aziz, supplicant, Baz, prayer space, omission, prostration, mosque, Abdullah", ## 5
                 "wash, one who fasts, fasting, fasting, to break fast, Ramadan, travel, dirty, Ramadan retreat, ritual cleansing", ## 6
                 "veil, youth, (sheikh) Tamim, Azzam, tanks, finery, wear, r(typo), sister, Peshewar", ## 7
                 "tithing, divorce, banks, divorce, card, banks, to perform pilgrimage, poor, sell, dirhams (currency)", ## 8
                 "Europe, civilization, European, mankind, church, goods, generations, their lives, huge, legend", ## 9
                 "to forbid, analogy, permission, general, evidence, forbid, text, absolutely, voting, concensus", ## 10
                 "excommunicate, apostate, apostasy, sponsorship, idolatry, excommunication, idols, to make permissible, our readers, apostasy", ## 11
                 "heterodoxy, innovator, Sufi, Salafi, to draw near to, distinguish, (the) saved (group), to undertake, innovation, multiply", ## 12
                 "Shari\'a, to legislate, to send down, to judge, judgment, justice, parliament, court, democracy, those", ## 13
                 "creatures, characteristic/quality, names, throne, al-Tahawi, proof, mean, question (abbreviation), heavens, explainer", ## 14
                 "narrate, Daoud, chain of narration, narrate, al-Bayhaqi, al-Tabrani, our saying, al-Tirmidhi, Father, (Abu) Hurayra") ## 15

## With all the words translated, I revised my topic labels
## I first printed out the translated words together, just for convenience
for(i in 1:15){
  print(paste("TOPIC",i))
  print(c(englishF[i],englishFREX[i]))
}

## I was still a bit puzzled for some of them: 2, 7, 9, 12, 14
## I used findThoughts to identify documents that were representative of each topic
## Below are my notes while reading through.

## topic 2
findThoughts(mod.out, topic=2, texts=as.character(meta$path),n=10)$docs
## the top documents are social theory by Muhammad and Sayyid Qutb.  I'm calling it Social theory

## topic 7
findThoughts(mod.out, topic=7, texts=as.character(meta$path),n=10)$docs
## although the FREX terms come from Abdullah Azzam on this one, the top documents are fatwas on family/women's issues

## topic 9
findThoughts(mod.out, topic=9, texts=as.character(meta$path),n=10)$docs
## similar to the Social theory top documents, but I'm going to call this "islam and modernity"

## topic 12
findThoughts(mod.out, topic=12, texts=as.character(meta$path),n=10)$docs
## the top documents are about Salafism -- seems to be about proper practice of Islam

## topic 14
findThoughts(mod.out, topic=14, texts=as.character(meta$path),n=20)$docs
## The top documents are all sections of a book on creed.  The fatwas are also about Creed -- the nature of god, etc.

## Getting exemplar documents for the Jihadist topics
## topic 11
findThoughts(mod.out, topic=11, texts=as.character(meta$path),n=20)$docs
findThoughts(mod.out, topic=1, texts=as.character(meta$path),n=30)$docs

## I'm declaring the topics labeled.  These are the final names.  It took me 6.5 hours or so to do this whole process.
topicnames <- c("Fighting","Social theory","Politics","The Prophet","Prayer",
                "Ramadan","Family and Women","Money, Pilgrimage,\nand Marriage","Islam and Modernity","Hadith",
                "Excommunication","Salafism","Shari\'a and Law","Creed","Hadith Narration")

################################################################
## make a regression table
################################################################

## simple regression plot
prep1 <- estimateEffect(1:15 ~ combinedjihadi, mod.out,meta)
prep <- plot.estimateEffect(prep1, "combinedjihadi", model=mod.out, method="difference",cov.value1=1,cov.value2=0,custom.labels=topicnames,labeltype="custom")


## make a custom regression plot (the one in the paper)

## save objects with the topic numbers
topics <- 1:15
revtopics <- rev(topics)
revtopics <- rev(seq(1,15*1,1))

## this is a custom graphing parameter I made for myself -- it's used below
offset <- 4.5
## set the working directory to the output directory
setwd(outdir)
## start the image file
png("arabresults.png",width=1000*7, height=1100*7, res=150*7, pointsize=8)
## set the graphing parameters
par(mar=c(3,2,3,1))
## get the confidence intervals
xlow <- unlist(lapply(prep$cis,function(x){return(x[1])}))
xhigh <- unlist(lapply(prep$cis,function(x){return(x[2])}))
## get the estimates
est <- unlist(lapply(prep$means,function(x){return(x[1])}))
## set the graph labels to blank
xlab <- ""
ylab <- ""
main <- ""
## start the plot
plot(0, 0, col = "white", yaxt = "n", xlab = xlab, main = main, axes=F,
         xlim = c(min(xlow) - offset, max(xhigh) + 0.15), ylim = c(.5,
                                                        max(revtopics)), ylab = ylab, pch=16)
## make alternating gray and white background
ss <- revtopics[which(1:length(revtopics) %% 2==1)]
gap <- (revtopics[1]-revtopics[2])/2
for(jj in ss){
      polygon(x=c(-100,100,100,-100),y=c(jj-gap,jj-gap,jj+gap,jj+gap), border=NA,col="gray90")
}
## spread out the estimates so that differences are visible
spread <- 2.5  ## larger values spread out the coefficient points
## add the axis for the coefficients
axis(1,at=c(-.5,0,.5),labels=c(-.5,0,.5)*spread)
axis(3,at=c(-.5,0,.5),labels=c(-.5,0,.5)*spread)
## make a line up the center
lines(c(0, 0), c(0, (max(revtopics) + 2)), lty=2)
## add the estimates
points(est*spread,y=revtopics, pch=20, cex=.75)
## add the confidence bands
segments(x0=xlow*spread,x1=xhigh*spread,y0=revtopics,y1=revtopics)
## add the topic #s and labels
text(rep(min(xlow),15) - (offset+.4), revtopics, paste(topics, ".\n", sep=""),cex=.75,font=2,pos=4, xpd=NA)
text(rep(min(xlow) - (offset+.2),15), revtopics, paste(topicnames, "\n", sep=""), cex=.75,font=2,pos=4, xpd=NA)
## add the words 
k <- 8  ## this is the number of words to use -- I use 8 because it fits nicely
## put the words into a holder
holder <- c()
for(i in 1:15){
  holder[i] <- paste("F:",paste(strsplit(englishF[i],", ")[[1]][1:k],collapse=", "),
    ## this line has a special k to make the words fit
  "\nFREX:",paste(strsplit(englishFREX[i],", ")[[1]][1:min(c(k,length(strsplit(englishFREX[i],", ")[[1]])))],collapse=", "),
  "\nF:",paste(topicslists$prob[i,1:k],collapse=", "),
  " FREX:",paste(topicslists$frex[i,1:k],collapse=", ") )
}
## Then add the words from the holder into the plot
text(rep(min(xlow) - (offset-.6),15), revtopics, holder, cex = 0.75, pos=4)
## finish plotting
dev.off()


#############################################################
## scatter plot of excommunication and fighing topics
#############################################################

## I use this library for "rename.vars()"
library(gdata)
## get out the topic proportions for the fighting topic
meta$stm.est.jihadi<-mod.out$theta[,fTopicNum]
temp<-aggregate(meta$stm.est.jihadi,by=meta["cleric"], mean, na.rm=TRUE)
temp<-rename.vars(temp,from="x",to="jihadi")
## get out the topic proportions for the excommunication topic
meta$stm.est.excom<-mod.out$theta[,eTopicNum]
temp2<-aggregate(meta$stm.est.excom,by=meta["cleric"], mean, na.rm=TRUE)
temp2<-rename.vars(temp2,from="x",to="excomm")
## make an object with the topic proportions for each author
topics.clerics<-cbind(temp,temp2[,2])
topics.clerics<-rename.vars(topics.clerics,from="temp2[, 2]",to="excomm")

## get the metadata
setwd(datadir)
dat <- read.csv("jihad_metadata_edited.csv",as.is=T)
dat <- dat[dat$combinedjihadi!=9,]
## set the rownames of "topics.clerics" to be the cleric names
rownames(topics.clerics)<- topics.clerics$cleric
## save it as a temporary object before combining
tmp <- topics.clerics[dat$cleric,2:3]
## set the column names
colnames(tmp) <- c("fighttopic","excommtopic")
## combine the metadata and the author topic proportions
dat <- cbind(dat,tmp)
## subset to just the variables we want
dat <- dat[,c("cleric","fighttopic","excommtopic","combinedjihadi","yearofbirth")]
## uniquify
dat <- unique(dat)

## make the scatter plot for the paper
setwd(outdir)
## start the image file
pdf("jihadtopics.pdf",7,5)
## plotting parameters
par(mar=c(5,4,1,1))
plot(dat$excommtopic,dat$fighttopic,type="n", xlab="Excommunication topic proportion",ylab="Fighting topic proportion",
  xlim=c(-.15,.45),ylim=c(-.01,.4),axes=F)
## add the axes
axis(1);axis(2,las=2)
## make a gray background
polygon(x=c(-100,100,100,-100),y=c(-100,-100,100,100),col="gray90")
## add the white lines
abline(v=seq(-.5,1,.1),col="white")
abline(v=0,lwd=3,col="white")
abline(h=seq(-.5,1,.1),col="white")
abline(h=0,lwd=3,col="white")
## add the data points
points(dat$excommtopic[dat$combinedjihadi==0],dat$fighttopic[dat$combinedjihadi==0],col="blue",pch=2)
points(dat$excommtopic[dat$combinedjihadi==1],dat$fighttopic[dat$combinedjihadi==1],col="red",pch=19)
## identify individual data points (commented out)
#identify(dat$excommtopic,dat$fighttopic,dat$cleric)
## Add the names of individual data points
text(dat$excommtopic[dat$cleric=="asam0.bn.ladn"],dat$fighttopic[dat$cleric=="asam0.bn.ladn"],"Usama bin Laden",pos=1,col="red")
cc <- "syd.QTb"
text(dat$excommtopic[dat$cleric==cc],dat$fighttopic[dat$cleric==cc],"Sayyid Qutb",pos=2,col="red")
cc <- "3ly.alKDyr"
text(dat$excommtopic[dat$cleric==cc],dat$fighttopic[dat$cleric==cc],"Ali al-Khudayr",pos=4,col="red")
cc <- "a7md.alKaldy" 
text(dat$excommtopic[dat$cleric==cc],dat$fighttopic[dat$cleric==cc],"Ahmed al-Khalidi",pos=2,col="red")
cc <- "3mr.3bdalr7mn"
text(dat$excommtopic[dat$cleric==cc],dat$fighttopic[dat$cleric==cc],"Umar Abd-al-Rahman\n(the blind sheikh)   ",pos=3,col="red")
cc <- "3bd.allh.3zam"
text(dat$excommtopic[dat$cleric==cc],dat$fighttopic[dat$cleric==cc],"Abdullah Azzam",pos=2,col="red")
cc <- "abw.ala3lA.almwdwdy"
text(dat$excommtopic[dat$cleric==cc],dat$fighttopic[dat$cleric==cc],"                     Abu al-Ala al-Mawdudi",pos=3,col="red")
## Add the legend
text(.28,.4,"Jihadists",xpd=NA,col="red",cex=1.2,pos=4)
text(.28,.37,"Non-Jihadists",xpd=NA,col="blue",cex=1.2,pos=4)
points(.28,.4,col="red",pch=19)
points(.28,.37,col="blue",pch=2)
## end the plot
dev.off()


#############################################################
## Topic network figure
#############################################################

## set the directory to the output dir
setwd(outdir)

## This is the matrix of Correlations Between Topics
cor(mod.out$theta)

## Size of Topic: Size depends on how you calculate it.  mod.out$theta is a D-by-K matrix with document d in 1:D and its loading onto each topic.  If you want frequency by word tokens then you just have to multiply through the word counts within each document.
head(mod.out$theta)
## I use this vector of wordcounts
wordcounts <- unlist(lapply(data, function(x) sum(x[2,])))
## there are fractional wordcounts due to variational approximation.
round(mod.out$theta[,1] * wordcounts,2)

## Calculate the proportion of words devoted to topics
topicPropsInCorpus <- rep(NA,15)
for(i in 1:15){
  topicPropsInCorpus[i] <- (sum(mod.out$theta[,i] * wordcounts))/sum(wordcounts)
}
## This now holds the topic proportions in the corpus
topicPropsInCorpus
## sums to one, as it should
sum(topicPropsInCorpus)
## add the topic labels
names(topicPropsInCorpus) <- topicnames

## Plot the network
library(igraph)
## using a non-binary distance matrix
mat2 <- cor(mod.out$theta)
## setting the negatives to zero
mat2[mat2<0] <- 0
## setting the diagonal to zero
diag(mat2) <- 0
## this gives us positive correlations between topics
mat2
## rename this object as "out"
out <- mat2

## set a seed so that it's reproducable
set.seed(2138)
## make the graph object
 g <- graph.adjacency(out, mode="undirected", weighted=T)
  if(length(labels)==0) labels = paste("Topic", topics)
## make the edges thickness proportional to correlation
cor(mod.out$theta)[,1]
E(g)
edges <- get.edgelist(g)
edgecors <- rep(NA,nrow(edges))
for(i in 1:nrow(edges)){
  edgecors[i] <- cor(mod.out$theta)[edges[i,1],edges[i,2]]
}
edge.width=35*edgecors
## look at and set other graph parameters
E(g)$weight
E(g)$size <- 1
E(g)$lty <- 1
E(g)$color <- "black"

#Make the colors indicate the direction of the coefficient
## get the estimates
est <- unlist(lapply(prep$means,function(x){return(x[1])}))
## get colors
mycols <- rev(colorRampPalette(c("red", "white", "blue"), bias=1)( 20 )) ## (n)
## assign the color category for each coeff
seq(-.3,.3,length.out=21)
colcat <- rep(NA,length(est))
for(i in 1:length(colcat)){
  colcat[i] <- max(which(est[i] > seq(-.3,.3,length.out=21)))
}
## These are now the color category for each coefficient
colcat
## and these are the associated colors
mycols[colcat]
## This checks to make sure the colors are working as I expect
## Blue should be on the left and red on the right.
plot(est,1:15,pch=19,cex=2,col=mycols[colcat]);abline(v=0)

## label the vertices
V(g)$label=topicnames
## set the size of vertices proportional to the proportion in the corpus
V(g)$size <- topicPropsInCorpus*200
## set the color of vertices
vertex.color = mycols[colcat]
## set other vertex characteristics
vertex.label.cex = 1
vertex.label.color = "black"
## set the edge color
edge.color = "gray60"
## set a seet so that the layout is reproduceable
set.seed(21)
## pull out the weights to include in the layout
wts <- E(g)$weight
## make the layout
mylayout <- layout.fruchterman.reingold(g,weight=wts)
## start the image file
pdf("corrNetwork.pdf",8,8)
## do the plot
plot(g, layout=mylayout,edge.color=edge.color,vertex.color=vertex.color, vertex.label.cex=vertex.label.cex, vertex.label.color=vertex.label.color,edge.width=edge.width)
## close the image file
dev.off()
## I then rearrange with photoshop to move the outlying vertices closer in


####################################################
## correlation of year and fighting + excommunication
####################################################

## calculate the proportion of words each author devotes to each topic
## make a holder
holder <- list()
for(i in 1:length(unique(meta$cleric))){
  holder[[i]] <- apply(apply(mod.out$theta,2,function(x){return(x*wordcounts)})[meta$cleric==unique(meta$cleric)[i],], 2, sum) / sum(wordcounts[meta$cleric==unique(meta$cleric)[i]])
}
names(holder) <- unique(meta$cleric)

## These are the topic profiles of some of the Jihadists who write much about fighting or excommunication
par(mar=c(5,10,2,2),mfrow=c(1,3))
plot(rev(holder[["abw.ala3lA.almwdwdy"]]), 1:15, axes=F); axis(1); axis(2,at=1:15,labels=rev(topicnames), las=2)
plot(rev(holder[["syd.QTb"]]), 1:15, axes=F); axis(1); axis(2,at=1:15,labels=rev(topicnames), las=2)
plot(rev(holder[["m7md.QTb"]]), 1:15, axes=F); axis(1); axis(2,at=1:15,labels=rev(topicnames), las=2)

## make the figure predicting topic usage as a function of year birth
## copy the metadata
dat2 <- dat[dat$combinedjihadi==1,]
## add a new variable to it that's the sum of the excommunication and fighting topics
dat2$bothJtopics <- unlist(lapply(holder,function(x){return(sum(c(x[fTopicNum],x[eTopicNum])))}))[dat2$cleric]
## start the image file
pdf("yearofbirth.pdf",6,4)
par(mar=c(5,6,1,1))
plot(dat2$yearofbirth[dat$combinedjihadi==1],dat2$bothJtopics[dat$combinedjihadi==1], xlim=c(1880,1980),ylim=c(0,.5), axes=F, 
  xlab="Year of cleric birth",ylab="Proportion of writing devoted to\nexcommunication and fighting")
axis(1);axis(2,at=seq(0,.5,.1),labels=T,las=2)
abline(lm(bothJtopics~yearofbirth,data=dat2))
## I used this to look at which clerics are which
#identify(dat$yearofbirth[dat$combinedjihadi==1],dat$bothJtopics[dat$combinedjihadi==1],dat$cleric[dat$combinedjihadi==1])
dev.off()


